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We investigate the evolutionary dynamics of an idealised model for the robust self-assembly of 
two-dimensional structures called polyominoes. The model includes rules that encode interactions 
between sets of square tiles that drive the self-assembly process. The relationship between the 
model's rule set and its resulting self-assembled structure can be viewed as a genotype-phenotype 
map and incorporated into a genetic algorithm. The rule sets evolve under selection for specified 
target structures. The corresponding, complex fitness landscape generates rich evolutionary dy- 
namics as a function of parameters such as the population size, search space size, mutation rate, 
and method of recombination. Furthermore, these systems are simple enough that in some cases 
the associated model genome space can be completely characterised, shedding light on how the 
evolutionary dynamics depends on the detailed structure of the fitness landscape. Finally, we apply 
the model to study the emergence of the preference for dihedral over cyclic symmetry observed for 
homomeric protein tetramers. 

PACS numbers: 61.46.Bc, 87.10.Mn, 87.23.Kg, 81.16.Dn 



I. INTRODUCTION 

Self-assembly processes, in which constituent compo- 
nents reliably assemble into a complete structure with- 
out external control, are ubiquitous in nature, providing 
the means by which sophisticated biological machinery 
such as protein complexes are formed within organisms 
PP. A key question then arises: how did the interactions 
that drive these self-assembly processes evolve over bil- 
lions of years to form the optimised systems we observe 
today [2M]? Bioinformatic studies of protein complexes 
[5] suggest that a number of observed trends in protein 
quaternary structure are caused not only by the biolog- 
ical function under selection, but also by the details of 
the evolutionary dynamics. Some of these trends have 
recently been explained by using computer simulations 
of a simple continuous patchy particle model [6] for glob- 
ular proteins [7] [S]- However, such models are compu- 
tationally expensive because a detailed simulation of the 
assembly process is required at each step in evolutionary 
time. 

In this paper we study the evolutionary dynamics of 
a highly idealised coarse-grained model for the evolu- 
tion of self-assembling systems, for which the assembly 
process can be simulated quickly and straightforwardly. 
The model consists of an 'alphabet' of square tiles that 
self-assemble into polyominoes: unions of connected cells 
on a 2D square lattice. The alphabet of available tiles, 
which we term the assembly rule set, contains a descrip- 
tion of the interactions that drive the assembling system 
towards a final structure [9 . A physical interpretation 
of the model consists of a structure assembling on a 2D 
substrate in contact with a suspension of tiles, as shown 
in Fig. [T] These tiles can form many kinds of structures, 




FIG. 1: (colour online) Illustration of a possible realisation of 
physical polyomino assembly. Square tile building blocks interact 
with each other through complementary bonding between edges, 
here illustrated with interacting polymer chains. In addition, tiles 
experience an attractive interaction to a flat substrate, leading to 
growing polyomino structures on a surface. 



both bounded and unbounded. We focus on determinis- 
tic rule sets that always assemble into the same bounded 
2D structures, a class of behaviour that is analagous to 
the monodisperse self-assembly observed for example for 
many kinds of protein quaternary structures. 

These models may also be relevant for experimental 
systems such as 2D self-assembled systems that have been 
made of RNA [10] and DNA [H] tiles. Each tile can be 
tailored to interact with its neighbours through comple- 
mentary bonding. Patterns and grids of varying geome- 
tries on the nanoscale have been produced by changing 
these design rules, with some examples being circuit pat- 
terns T2] and Sierpinski triangles [T3]. The variety of 
structures that can be produced using DNA tiling [T3] 
and DNA-linkcd particles [15] is rapidly increasing. The 
evolutionary design of polyomino structures may shed 
light on the design of these synthetic systems. 
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Pioneering work by Wang [T5] demonstrated that tiles 
could be used to specify a Turing machine. In an impor- 
tant development, Winfree showed that DNA nanotech- 
nology could be used to create molecular Wang tiles [TT] . 
Self-assembling tile sets can thus perform computational 
tasks such as binary counting, and a measure of the com- 
plexity of assembly sets required for such algorithmic 
applications have been computed |17j . This theoretical 
work has been extended to study the details of tile as- 
sembly nucleation [18] and the effects of errors in the 
assembly process [ft)] . 

In this study, we use genetic algorithms (GAs) [2"0]l2"T] 
that search through the space of all possible rule sets 
to find those that generate the deterministic assembly 
of desired polyomino structures. Despite its simplicity, 
and resulting computational tractability, the model pro- 
duces rich evolutionary behaviour. The assembly process 
can be viewed as a mapping that transforms an assem- 
bly rule set into an assembled polyomino structure. This 
mapping is reminiscent of the genotype-phenotype map in 
evolutionary biology, whereby information in the genome 
(the genotype) is used to develop the physical form of a 
biological structure (the phenotype). 

We investigate how the evolutionary dynamics of our 
model system depends on parameters such as popula- 
tion size, mutation rate and recombination. In GAs, mu- 
tation rate has been shown to dramatically affect the 
speed of evolution, with populations evolving at higher 
rates around an optimal mutation rate that is roughly 
the reciprocal of the genome length [531 H3] • Biological 
organisms also often have mutation rates around this op- 
timal value [24J [25] . Recombination has also been shown 
to increase the speed of evolution on a simple fitness 
landscape [26 . We study how these evolutionary vari- 
ables affect the adaptation and discovery times of our 
self-assembling systems. 

An important property of our model is that it is 
simple enough to allow, in some cases, an exhaustive 
search of the associated search space, yielding a fully- 
characterisable but highly non-trivial fitness landscape 
[27]. [28] that facilitates a detailed analysis of the under- 
lying evolutionary dynamics. 

In addition, we aim to explore the emergence of sym- 
metry in evolving self-assembling systems. It has been 
observed, for example, that homomeric tetramer protein 
complexes show a strong preference for dihedral (-D2) 
symmetry over cyclic (C4) symmetry [7]. We study 
this preference as a function of various evolutionary pa- 
rameters with our simplified polyomino system, for which 
a complete characterisation of the fitness landscape can 
be achieved. 

This paper is structured as follows. In Section|ll]we de- 
scribe our model of self-assembling polyominoes, and our 
implementation of genetic algorithms. In Section |III| we 
exhaustively study the search space defined by a partic- 
ular parameterisation of our model. Section |IV| analyses 
how evolutionary variables including mutation rate, pop- 
ulation size and search space size affect the dynamics of 
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FIG. 2: (colour online) Illustration of polyomino assembly, for a 
rule set with nt = 4, n c = 7, and with the nucleus and interaction 
conventions described in the text. A binary representation of a rule 
set is translated to nucleus, tile and interaction information. The 
nuclei are placed on a lattice in the initiation stage, and growth 
progresses stochastically, to a final output possessing shape and 
tile, but not orientational, determinism: the diagonal neighbours 
of the central tile have two possible orientations, as their 4 edge 
can bond to either of the two adjacent 3 edges. 



polyomino evolution. In Section [V] we apply our model 
to study the evolution of homomeric tetramer protein 
complexes, and we list our conclusions in Section |Vl] 



II. MODEL & METHODS 
A. Model Implementation 

Our model uses interacting square tiles to model the 
self-assembly of 2D polyomino structures on a square 
lattice [9 . The interactions between adjacent tiles are 
defined by the nature of each tile's edges, which are as- 
signed 'colours', with any two colours either experiencing 
no interaction or an attraction. In this conceptual model, 
there is no energy or temperature scale, so two edges are 
either non-interacting or have an effectively infinite at- 
tractive interaction, making bonding irreversible. 

A given assembly scenario will consist of n t tile types, 
and an alphabet of n c available colours. Each tile is en- 
tirely specified by a description of its four edge colours. 
We will denote a tile as an ordered set of four colours, 
with the first element corresponding to the top edge and 
subsequent elements corresponding to the edges reached 
in clockwise order, for example, {1,2,3,4}. An n c x n c 
binary interaction matrix A describes the interaction 
between colours, with colours i and j experiencing an 
attractive interaction if Aij = 1, and no attraction 
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FIG. 3: (colour online) (a) Unbound and (b) non-deterministic 
rule sets. The interactions in this case follow the convention in 
Eqn. [I] so edge types 1 and 2 experience an attractive interaction 
and is neutral. The rule set in (a) consists of one block that 
is capable of bonding to itself, thus creating an endless chain of 
repeated blocks. The rule set in (b) contains a block capable of 
bonding to more than one other, leading to tile non-determinism 
depending on which bonding block is added first. In this case, the 
two different bonding tiles create different structures, leading to 
shape non-determinism. 



studies, seed particles constructed of DNA form the nu- 
cleus of a structure and contain information to regulate 
the assembly process. This approach effectively corre- 
sponds to an SFN setup. 

We will adopt conventions for the nucleus tiles and 
the structure of the interaction matrix A, allowing us to 
simplify the representation of a rule set. We will use 
an SFN, and take the nucleus tile to be of the tile type 
first described in the rule set. Furthermore, we fix the 
orientation of the nucleus tile, so that the edge specified 
first in the rule set is taken to be the upper edge of the tile 
when first placed on the grid. Under our convention, the 
position of the nucleus tile is arbitrary, and polyominoes 
that differ only by translations are counted as equivalent. 

We will usually (with an exception in Section Vj which 
allows the incorporation of self-interacting colours) fix 
the interaction matrix by defining the interaction be- 
tween colours i and j (represented by non-negative in- 
tegers) as: 



otherwise. The generalised case of varying interaction 
strengths has been studied analytically [29] , but for sim- 
plicity we consider binary interactions. 

The tiles are similar to Wang tiles [311], with two im- 
portant differences: interactions between colours are not 
limited to each colour bonding only with itself, and the 
tiles may be rotated to any of the four possible orienta- 
tions allowed by C4 symmetry (for example, {1, 2, 3, 4} = 
{3, 4, 1, 2}). The sides of a tile therefore comprise what is 
termed 'an rt c -ary fixed necklace' of length 4 [9], [31] . The 
generalisation to free necklaces [31 , in which tiles may 
also be 'flipped' ({1,2,3,4} = {1,4,3,2} and its cyclic 
variants), will be visited in Section |v| 

Assembly progresses on an infinite square lattice, and 
takes places in two phases: initiation and growth (see 
Fig. [2| . The initiation phase involves one or more tiles 
being placed on the (initially empty) lattice at prescribed 
positions and orientations: these are the nucleus tiles, 
each of which is described by the tile type of the nucleus, 
its co-ordinates on the lattice, and its orientation. The 
combined instruction set representing nucleus, tile edge, 
and interaction matrix data is the rule set for a particular 
assembly scenario. 

There are several alternative schemes for nucleating 
assembly in this model. Assembly may progress from a 
single initial tile, laid down at the start of the assembly 
process. In this case, the single tile may be of a fixed, 
specific tile type — which we will term a single fixed nu- 
cleus (SFN) — or of a tile type arbitrarily chosen from 
the rule set — which we will term a single general nucleus 
(SGN). It has been shown that to guarantee determin- 
istic assembly from an arbitrary nucleus tile, consider- 
ably more information content is often required within 
genomes [9]. 

The question of nucleating tile-based self-assembly has 
been addressed theoretically [32] and experimentally [TH] 
in the context of algorithmic DNA assembly. In these 
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so that each colour only interacts with one partner, 1 O 
2, 3 o 4,... and provides a neutral edge, which does 
not interact with any other edge type. 

The combination of conventions for assembly nucle- 
ation (SFN, with the first tile specified in the rule set as 
the nucleus) and the interaction matrix (Eqn. [I]) allows 
us to represent a given rule set by specifying the edges of 
the tiles involved in assembly alone. Rule sets can then 
be represented straightforwardly by a binary string (see 
Fig. El), by writing each numerical parameter in the rule 
set (each tile edge) as its binary counterpart and con- 
catenating all the binary variables into one long string. 
This resulting 'genome' is then suitable for processing 
with genetic algorithms (see Section II C I . 

Growth progresses stochastically in the following man- 
ner. A tile type is chosen randomly from a uniform dis- 
tribution over the available tiles. A position on the lat- 
tice is selected randomly with the constraint that it must 
be adjacent to a previously laid tile. The chosen tile is 
cycled in random order through its four possible orien- 
tations at the chosen point. If during this cycling the 
tile experiences an attractive interaction to any of its 
four neighbouring lattice points, it bonds immediately 
in that configuration at the chosen site, as illustrated in 
Fig. [2] In this way, bonding occurs irreversibly, but the 
model can be generalised to allow reversible interactions 
by introducing a temperature scale, relaxing the binary 
constraint on interaction matrix A, and allowing assem- 
bly to proceed within a simulation that includes thermal 
effects. 



B. Classes of Assembly Behaviour 

Rule sets in our model may result in unbound struc- 
tures: those where the assembly process proceeds in at 
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FIG. 4: (colour online) Illustration of some UND polyomino types 
resulting from growth of genomes with nt = 2, n c = 8. UND 
polyominoes form the majority of achievable structures. The two 
colours label the two different tile types that may be involved in 
assembly. Letters in brackets denote whether the structures are 
bound (b), unbound (u), deterministic (d), non-deterministic (nd), 
space-filling (sf) and periodic in one (ID) or two (2D) dimensions. 



least one direction without termination. Unbound struc- 
tures may result, for example, from a set of one or more 
tiles that bonds to itself repeatedly, forming an endless 
chain of repeated units, as illustrated in Fig. |3](a). 

Self-assembly in biology may also yield unbound struc- 
tures. Proteinaceous structures consisting of extended 
sets of repeated units include helical protein filaments 
such as microtubules [33], actin filaments [34] and to- 
bacco mosaic virus [35] ; two-dimensional arrays such as 
S-layers [35] and purple membranes [37] ; and even three- 
dimensional crystals [38] , although some biological mech- 
anism must usually be present to regulate the size of these 
assemblies and prevent them being truly unbound [55] . 

Another assembly feature that may result from our 
model is non- determinism, whereby the same set of rules 
may lead to different structures forming in the growth 
phase. This non-determinism is due to the inherent 
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FIG. 5: 61.1% of the polyominoes that can be grown from genomes 
with nt = 2,n c = 8 are UND structures (X). The rest are bound, 
deterministic polyominoes. The number of genomes that encode 
for each polyomino are given - for sets of structures with identical 
values, each structure occurs the given number of times in genome 
space. Some structures are given names for ease of reference in the 
text. 



stochasticity in the assembly process. Non-determinism 
may arise when a tile edge is capable of bonding to more 
than one other tile edge, which may occur, for example, 
when the partner to a given edge type appears on more 
than one tile within the rule set, as in Fig. |3jb). 

Non-determinism may occur in several ways. Shape 
non- determinism is the least subtle form, whereby the 
overall shape of the produced structure (disregarding any 
detail of tile types, sides and orientations) differs stochas- 
tically in different assembly runs (see Fig. [3] (b) for exam- 
ple). Tile non- determinism occurs when the same over- 
all structure is achieved for all runs, but sites within the 
structure are occupied by different tile types stochasti- 
cally. Orientational non- determinism occurs when the 
structure is both shape- and tile-deterministic, but tiles 
within the structure differing stochastically in orientation 
between assembly runs (an example of this is the struc- 
ture in Fig. [2]). Another type of non-determinism, steric 
non- determinism, may also occur as a result of the differ- 
ent speeds of growth in two directions that converge on 
the same point: if two arms of a structure pass through 
the same lattice point, the structure will differ depend- 
ing on which arm arrives there first and hinders growth 
of the other. This type of non-determinism does not re- 
quire the multiple bonding edges mentioned above, and 
is thus hard to detect through observation of the genome. 

In biology, non-determinism can also occur in a num- 
ber of ways. Some closely related proteins coassemble 
into complexes of well-defined size and shape, but in 
which the identity of the protein at any position is ran- 
dom. An example of this phenomenon is in the seeds 
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of pea plants [30] where the legumin protein is formed 
by a number of paralogous genes, which result in hex- 
amers containing randomly assorted subunits of similar 
but distinct polypeptide sequences. This example would 
correspond to tile non-determinism in our model. There 
also exist examples of shape non-determinism in biology, 
where proteins, such as certain heat shock proteins, as- 
semble into clusters with a polydisperse distribution of 
sizes gH 133. 

Finally, our assembly model may yield structures that 
are bound (of finite size) and deterministic (in which the 
self-assembly process always forms the same structure, 
with a specific shape). The majority of protein quater- 
nary structures fall into this category [55] . 

For completeness, we note that there is incomplete 
overlap between the sets of non-deterministic and un- 
bound structures: rule sets may code for outputs 
that are unbound but deterministic, bound but non- 
deterministic, bound and deterministic or unbound and 
non-deterministic. In this study, we will focus on bound, 
deterministic structures, and will refer to structures not 
meeting these criteria as UND structures (unbound or 
non-deterministic). Some examples of these structures 
are shown in Fig. |4j 

All these forms of non-determinism can in theory be 
detected by running each growth phase a large number 
k times and comparing the output each time. We shall 
employ k = 10, a value that was confirmed through pre- 
liminary investigation to detect most non-deterministic 
structures while retaining computational speed. In this 
investigation, we choose tile and orientational determin- 
ism as our desirable criterion. 



C. Genetic Algorithm Details 



individuals preserved, immune to the effects of mutation 
|21) . We will explore the use of elitism in Section IV F 
but will generally set e = 0. 

GAs may employ crossover, modelling recombination. 
Without crossover, in the asexual regime, new individ- 
uals begin as cloned copies of selected genomes. With 
crossover, modelling sexual reproduction, new individu- 
als are formed by selecting two 'parent' rule sets from 
the old generation, forming a new rule set by combining 
the rule sets of these parents. The crossover scheme we 
employ is single-point crossover, where the first Lr bits 
from one parent and the last L — Lr bits from the other 
are combined to form a new individual, and Lr is chosen 
randomly from [0, L]. 

The implementation of crossover in a simulation is con- 
trolled by the crossover rate R, giving the proportion of 
new genomes that are formed through crossover. For sim- 
plicity, we will only employ values of R = (correspond- 
ing to asexual reproduction) and R = 1 (corresponding 
to sexual reproduction). 

Another genetic operator used in GAs is mutation, 
whereby individuals in a generation undergo stochastic 
changes to their rule sets. We employ point mutation, 
whereby each bit in the genome is flipped with probabil- 
ity /!. 

Genomes may contain redundant information, with a 
tile type being coded for more than once in the binary 
string. In addition, information on tiles and edges that do 
not play a role in the assembly of the final structure may 
be included in the genome. This unused information in 
genomes allows neutral mutation to progress. A genome 
may also, in the aforementioned non-deterministic case, 
code for many different polyomino structures, and the 
same structure may be produced by more than one 
genome, providing a many-to-many mapping. 



Genetic algorithms (GAs) are a class of optimisation 
procedures that employ operators based on evolutionary 
biology to reach a solution to some problem I2"T] . 
Typically, GAs involve a population of N individuals, 
each representing a trial solution to a problem. A fitness 
function quantitatively measures the performance of an 
individual at solving the required problem. 

GAs take place over a number of time steps or gen- 
erations. Each generation, the fitness function is used 
to assign a fitness value fa to each genome j in a pop- 
ulation of N individuals. A selection operator is then 
applied, selecting genomes for reproduction according to 
their fitness values, with high-fitness genomes being pref- 
erentially selected. The N rule sets comprising the next 
generation are then formed from selected genomes. We 
employ the roulette- wheel selection method [43], where 
the probability P(i) of a genome i being selected is pro- 
portional to its fitness: P(i) = fi/Ylj fj- 

A common practise in the implementation of GAs is 
to preserve a certain number of the fittest individuals 
in a population from one generation to the next. This 
approach is termed elitism, with a proportion e of fit 



III. SEARCH SPACE ANALYSIS 

The process of evolution can be viewed as an optimisa- 
tion process on the high-dimensional search space defined 
by all possible genomes [HZIEH] - An advantage of our self- 
assembly model is that the search space for simple struc- 
tures can be fully characterised. We first investigate the 
structure of the search space for a polyomino model with 
two tiles (n t = 2) and up to eight colours (n c = 8), allow- 
ing three bond types (1o2,3o4,5h6), with colours 
and 7 corresponding to neutral edge types. Each of the 
8 tile edges can be represented by log 2 8 = 3 bits, giving 
a binary genome of length L = 24. The search space 
therefore consists of 2 24 ~ 1.6 x 10 7 individuals. We will 
refer to search spaces as 5 rat ,n c , labelled by the number 
of blocks (tiles) n t and number of colours n c , so that the 
aforementioned search space is 52,8- 

We adopt the convention that the first tile encoded 
in the genome is the assembly nucleus, and its initial 
orientation is specified by the order in which its edges 
are encoded, with the top edge first and others following 



6 




FIG. 6: (colour online) Transitions between polyominoes in the 52,8 [fH = 2, n c = 8) system, (a) The value of a pixel denotes the total 
number of single-point mutations that result in a change from phenotype x to phenotype y, over all genotypes in 52,8 that encode x and 
y. (b) The value of a pixel denotes the average proportion of mutations that cause an x — > y transition, where the average is taken over 
all occurrences of x in 52,8- Pink pixels denote transitions between phenotypes that cannot be accomplished with a single mutation. 



in a clockwise direction. We then exhaustively evalu- 
late all polyomino structures that may be constructed 
in this system. The majority are UND structures, some 
examples of which are shown in Fig. [4] to illustrate the 
diversity of achievable forms. These structures include 
non-deterministic, bound structures (for example, A in 
Fig. [2]) , deterministic structures that are translationally 
periodic in one (D) or two dimensions (E, K), some of 
which may be space- filling (G). Unbound structures dis- 
playing shape- but not tile-determinism order also exist 
(F, M). 

The resulting structures are illustrated, along with the 
volume of search space they occupy, in Fig. [5] Sets of 
genomes encoding the same phenotype form the neutral 
network of a given phenotype. The large differences in 
neutral network size corresponding to different pheno- 
types are related to the differing amounts of information 
required to specify bonds for different phenotypes. For 
example, the single tile phenotype only requires an ab- 
sence of any bonding edges rather than any specific inter- 
action pairs, and correspondingly occupies a large volume 
of genome space. By contrast, the 4x4 block phenotype 
requires two interacting pairs of edges, at specific posi- 
tions relative to each other, and the number of genomes 
fulfilling these criteria is much smaller. 

In addition, all single mutation transitions were 
recorded, identifying the effect of every possible single 
mutation on every possible genome — which may change 
the phenotype or be neutral (with no phenotypic effect). 
Fig. |6|a) depicts the number of possible single-mutation 
transitions between different phenotypes, whilst Fig.J^b) 
depicts the probability of a transition to another struc- 
ture given an initial structure. In (b), the total number 
of transitions between two phenotypes are normalised by 
the number of genomes encoding the x-axis phenotype 



(see Fig. [1]). The resulting quantity measures the aver- 
age number of mutations in a genome that encodes phe- 
notype x that cause a transition to phenotype y. 

The Fiedler eigenvalue method [H] was used to arrange 
the phenotypes in Fig. [6] to maximise the "blockiness" 
of the resulting matrix by clustering rows and columns 
whose elements follow similar trends. This method no- 
ticeably groups modularly-related polyominoes — for ex- 
ample, the 2x1 and 3x1 structures are clustered to- 
gether, as several single-mutation changes allow transi- 
tions between these structures through the addition or 
subtraction of a single block. This clustering reflects the 
fact that pairs of polyominoes that share modules (tiles 
or bond sequences) are more closely connected in genome 
space than unrelated structures. 

Fig. [6^b) shows that the majority of single mutations 
from a given phenotype are either neutral, preserving the 
phenotype - leading to high diagonal values in the plot - 
or cause a transition to a UND or single-tile phenotype. 
The fraction of neutral mutations is noticeably smaller 
for larger polyominoes (for example, the 'Catherine wheel' 
structures and the 4x4 block have diagonal values under 
0.3) than smaller ones (for example, the single tile, 2x1 
blocks and the 2x2 block have diagonal values over 0.6), 
partially because genomes encoding small polyominoes 
contain more redundant information than those encoding 
large polyominoes. 

Another observable feature of the search space is that, 
for several phenotypes, the most common result of mu- 
tations that are not neutral and do not result in a UND 
phenotype is a loss of part of the structure associated 
with the phenotype. For example, a significant propor- 
tion of mutations lead from the 4x4 block to the 2x2 
block, removing the outer 'shell' of tiles. The T-shaped 
tetrominoes also show many transitions to the L-shaped 
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FIG. 7: (colour online) Fitness curves during a typical evolution- 
ary run. A population of genomes is evolved towards a structure 
of size s > 16 using the fitness function in Eqn. [2] at fiL = 0.5, 
R = 0, N = 10. The plot shows the mean (dashed) and maximal 
(solid) fitness within a population as time progresses. 



triominoes, as one tile is lost from the phenotype. These 
triominoes in turn show many transitions to the 2x1 
blocks, from the loss of another bonded tile. 

Fig. [6jb) gives a measure of the average robustness 
and evolvability [45] of a given phenotype. The diagonal 
values give the phenotypic robustness, measuring the av- 
erage (over all genomes that encode a given phenotype) 
number of possible mutations that preserve phenotype. 
This averaging gives a mean phenotype robustness rather 
than the robustness value for any individual genome |46) . 
Phenotypic evolvability can be measured in two different 
ways. Firstly, a sum over off-diagonal values gives the 
number of mutations that result in a useful (non-UND) 
phenotypic change. Secondly, the number of non-zero 
off-diagonal values in a column give the number of dif- 
ferent phenotypes that can be accessed from the source 
phenotype. The first measure can be used to describe 
the probability that a non-neutral mutation will result 
in a useful phenotype. The second is more closely re- 
lated to Wagner's definition of phenotype evolvability 
|46j : it measures the diversity of phenotypes accessible 
from the neutral network of a given phenotype. In our 
model, robustness and evolvability are related differently 
in different phenotypes: the Catherine wheel structures 
are highly evolvable according to both the above defini- 
tions, but have low robustness (about 0.3), whereas the 
2x2 square has high evolvability and high robustness 
(about 0.6). 



IV. EVOLUTIONARY DYNAMICS 

A. Evolving Polyomino Size 

In evolution, selection drives a system towards high- 
fitness phenotypes (analogous to a thermodynamic drive 



towards low-energy structures), and entropic effects 
favour those structures that occupy a large proportion 
of search space. This interplay of fitness and entropic 
terms is analogous to the concept of free energy in ther- 
modynamics, and indeed several studies have analysed 
evolution using a 'free fitness' quantity [47l|48j. It may 
be expected that the importance of a given phenotype 
in evolutionary dynamics is related to several factors, in- 
cluding the fitness of the phenotype and how frequently 
genomes that produce it occur in the search space. For 
example, if fitness is defined as proportional to polyomino 
size, we may expect large structures that occupy a large 
volume of search space (i.e. with relatively large neutral 
networks) — like the 'Catherine wheels' and the 4x4 block 
in Fig. [5] — to play important roles in evolutionary path- 
ways. 

Having characterised the 52,8 search space in detail, 
we now proceed to simulate evolution on a fitness land- 
scape in this search space, with a particular aim being to 
relate the evolutionary dynamics back to the structure 
of the underlying search space. We use a specific fitness 
function to drive evolution towards a given target with 
a GA. A simple example is evolution towards a bound, 
deterministic polyomino matching or exceeding a certain 
size, using the fitness function: 



1. 



F(P 1 ,P 2 ,...,P k ;s*) 



s(-Pi) > s*and all Pi 
identical and bound; 
s(P 1 )/s*, 

s(Pi) < s*and all P 
identical; 



0. 



Pi unbound orPi ^ Pj 
for any i,j. 



(2) 

Here the fitness function takes a set of polyominoes 
{Pi, ...,P C } produced through k repeats of the assembly 
process, and a desired size s* . The function returns a 
zero fitness value if the set of polyominoes is UND, and 
a fitness value proportional to polyomino size for bound, 
deterministic structures. A value of one means that a 
solution matching the size criterion has been found. 

We note that the previous section suggests that only 
a small minority of the possible mutations to a genotype 
lead to a phenotype of larger size. However, it may be ex- 
pected that on the rare occasions that such mutations do 
take place, selection will allow these phenotypic changes 
to be retained and propagate through the population. 

Fig.[7]shows the time evolution of a population of poly- 
ominoes towards the target s* = 16. On the £2,8 land- 
scape, only one phenotype fulfills this criterion: the 4x4 
block. We employ what we will term zero initial condi- 
tions, in which every bit in every genome at the start of 
the simulation is set to zero. In the self-assembly imple- 
mentation described above, this approach means every 
initial genome encodes a single tile phenotype, which is 
laid down and incapable of further bonding. 
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FIG. 8: (colour online) Adaptation time r (solid lines) and dis- 
covery time tjj (dashed lines) in generations, in 52,8 evolving to 
s* > 16, with mutation rate fi, at different N and with R = 0. 



The simulation begins with the trivial, single-tile phe- 
notype, then quickly 'discovers' beneficial interactions, 
increasing the size of the largest phenotype in the pop- 
ulation hrst to two then to four, with the 2x2 square 
structure being discovered. The mean fitness lags be- 
hind the maximal fitness, as many members of the pop- 
ulation will still possess lower fitness values - the mean 
fitness rises only gradually above the value correspond- 
ing to the single tile phenotype, as the information for 
the 2x2 square structure does not immediately prop- 
agate through the whole population. The slow spread 
of information is due to both the finite fitness advan- 
tage resulting from the larger size of the square structure, 
and the possibility of further mutations leading to UND 
structures. After several generations, a further beneficial 
interaction is discovered, creating the 'Catherine wheel' 
octomino, and in the next generation this structure is 
expanded upon to form the 4x4 structure. Note that 
the Catherine wheel structure is one of only a few phcno- 
types exhibiting a single-mutation transition to the 4x4 
block (see Fig. [6]). 

The discovery of the 4x4 block leads to a sharp rise in 
the mean fitness, which lasts several generations before 
flattening. This flattening is due to the non-zero muta- 
tion rate and the high transition probability between the 
4x4 block and other structures of lower fitness. This 
mutational entropy means that for a finite /z, perfect 
adaptation is not reached in this system. The simula- 
tion is terminated when more than half the population 
has maximal fitness. 



B. Varying Mutation Rate 

In GA experiments, the discovery time To measures 
how long a system takes to produce a single copy of a 
maximally fit solution, giving an indication of the speed 
at which evolution progresses. Specifically, To is the first 
generation in which at least one genome encoding a max- 
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FIG. 9: (colour online) Fitness curves with time for simulations in 
S 2 ,8 evolving to s* > 16, with N = 80, R = and (a) = 0.1, (b) 
fiL = 0.5, (c) i±L = 2. Maximum fitness in the population is shown 
in red (upper curves) and mean fitness is shown in green (lower 
curves). Adaptation, defined as the point where 50% or more of 
the population has maximal fitness, occurred at generations 95 for 
(a) and 51 for (b). (c) failed to adapt within the 20 000 generation 
cutoff. 



imally fit solution is present. 

The distribution of to in an ensemble of GA experi- 
ments is generally observed to be long-tailed, with infre- 
quent occurrences of very high discovery times. Due to 
computational limitations, we generally employ a cutoff 
of 20 000 generations in our GA runs. As these rare, high 
values can skew the mean of such a distribution, we use 
the median of the distribution as a measure for to , as this 
statistic is less prone to skew from the rare events and 
artefacts from the imposed cutoff. 1 000 GA runs were 
performed for each data point in the following plots. 

We measured the value of To in the 52,8 system, as a 
function of mutation rate \x at a range of population sizes 
N. We set R = and use zero initial conditions. Fig. [8] 
shows the results, tjj decreases monotonically with /i 
except in the case of low N, where a slight increase at 
high n is observed. The decrease in to at high fi is due to 
the allowed larger steps across search space and a more 
explorative search. The slight increase in tjj at high // in 
the low N case may be due to the inability of completely 
random search to efficiently explore the search space with 
a small population - in other words, either some memory 
of previously discovered information or a large population 
is required for optimal search. We will see in Sect ion [TV C| 
that the monotonic decrease in to for larger population 
sizes is due to the small size of the £2,8 search space, 
and that To exhibits an optimum with \x in larger search 
spaces. 
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Another timescale of interest in evolutionary simula- 
tions is the adaptation time r of a system, measuring 
how long a solution, designated as maximally fit, takes 
to dominate the population. We measure this quantity 
as the first generation in which more than half the popu- 
lation has maximal fitness. The reason this criterion was 
chosen is that, due to the high proportion of deleteri- 



ous mutations that decrease fitness (see Section III I , full 



adaptation is unlikely to occur in reasonably sized pop- 
ulations at finite \i due to the likelihood of at least one 
phenotype-changing mutation occurring in a population. 

Fig. [8] shows r values for the £2,8 system, with R = 0. 
A general observation is the presence of an optimal mu- 
tation rate ji* , at which t is a minimum. The optimal 
mutation rate arises from the following competition. At 
very low /x, r increases divergently as fi decreases. This 
increase in r at low fi is steeper at low N than at high 
N. The reason is simply that at low mutation rates, it 
takes a long time for the system to discover new phe- 
notypes, and this is made worse in smaller populations. 
On the other hand, arguments from population genetics 
[4"9"| suggest that full adaptation of a population becomes 
increasingly difficult for fi > 1/L, due to mutational en- 
tropy. Thus one expects an optimal mutation rate for 
adaptation around fi m 1/L. 

Fig. [9] shows examples of the time evolution of the 
fitness during simulations at a range of [i values (low: 
[iL = 0.1, intermediate: [iL = 0.5, high: fiL = 2). At 
low fi, the mean fitness closely tracks the maximal fit- 
ness, as diversity is low and the population is confined 
around a small region of genome space. The behaviour is 
due to the high correlation between generations: as little 
change is introduced to the gene pool through mutation, 
diversity in the population is low. 

At high /i, the mean fitness fluctuates around a low 
value, dominated by the entropic drive towards common, 
low-fitness structures (as most mutations are deleterious 
- see Fig. [6]). In this regime, the population is decor- 
related, and highly genetically diverse — resembling a 
random search across genome space. 

Behaviour at \xL = 1 is intermediate between these 
regimes, with some diversity resulting in a rather lower 
mean fitness than maximal fitness, but a clear relation- 
ship between the two showing that information is not 
being lost through decorrelation. 

The relationship between mean and maximal fitness 
also depends on the robustness of the phenotypes within 
a population. In Fig. [9] (b) , the mean and maximal fit- 
ness values are closer in magnitude for a local optimum 
(around generation 40) than for the global optimum (gen- 
eration 43 onwards). This difference suggests that the 
robustness of the global optimum is lower than that of 
the local optimum, as the population has more difficulty 
adapting to the fitter phenotype. 

We will use the terms exploration and exploitation to 
refer to the two regimes observable at high and low \i, 
respectively. Exploration refers to the random search 
regime at high /i, where genome space is explored uni- 



formly and randomly, and the entropic effect of muta- 
tion is too high for the population to become localised 
and adapt. Exploitation refers to the highly-correlated 
regime at low /1, where evolution progresses through 
small changes made to existing information, resembling 
a "hill-climbing" process with a low diversity. The inter- 
mediate fj, regime may be thought of as providing a com- 
bination of these two effects, with enough exploration to 
allow escape from local optima and enough exploitation 
to experience a drive to higher fitness values. 



C. Comparing Search Spaces 

To investigate the effect of changing the search space 
for the system, we next considered the Sq^ space, involv- 
ing rit = 6 blocks rather than the rit = 2 used previously. 
Genome length is now L = 72, with 2 72 ~ 4.7 x 10 21 
points in search space, more than 14 orders of magni- 
tude larger than the £2,8 space. We used a sampling 
approach, investigating 10 s points in Sq^, to investigate 
how the structure of this new search space may affect the 
search for an s > 16 structure. Firstly, a larger number 
of genomes in the new space encode for such a structure, 
with many possible ways of achieving the 4x4 square 
and other, more diverse structures with s > 16. How- 
ever, the associated exponential increase in the overall 
size of the search space means that a smaller proportion 
of genomes encode a structures with s > 16, with many 
more genomes now producing small or UND polyomi- 
noes. 

Fig. 10 shows the r and tjj behaviour with [i in 6>6.8- In 



this plot, we see first of all that even though the search 
space is many orders of magnitude larger, the optimal 
adaptation and discovery times are at most an order of 
magnitude larger. The qualitative behaviour of the dis- 
covery time td also shows an important difference from 
the simpler £2,8 system. This measure now exhibits an 
optimum with /z, generally around fiL > 1. At higher fi 
values, td increases, indicating that the large steps per- 
formed by high-/i search in this case are not beneficial. 
This optimum arises from a tradeoff between exploration 
and exploitation: the system must have a high enough 
/j to successfully explore a range of genome space, but 
must have a low enough fi so that useful information is 
not lost. 

At high fi, the gene pool decorrelates significantly from 
generation to generation, resulting in loss of information 
about intermediate-fitness structures that have been dis- 
covered. In the smaller £2,8 system, Fig. [8] suggests that 
this loss of information is not an important effect, as 
the highly random search afforded by high [i has a finite 
chance of discovering a suitable solution through explo- 
ration alone. However, in the exponentially-larger 56,8 
space, random search has a very low probability of dis- 
covering a suitable solution, and exploitation of existing 
information is important in the discovery of better solu- 
tions. 
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FIG. 10: (colour online) Adaptation time r (solid lines) and 
discovery time td (dashed lines) in 5s, 8 evolving to s* > 16, with 
mutation rate fi, at different N and with R = 0. 
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FIG. 11: (colour online) Adaptation time r (solid lines) and 
discovery time tx> (dashed lines) with random initial conditions in 
(a) 52,8 ancl (b) >56,8i evolving to s* > 16, with mutation rate fi, 
at different N and with _R = 0. 



D. Initial Conditions 



relevance when bits in a genome represent more funda- 
mental units of genetic information, as it corresponds to 
an interbreeding population with entirely different, ran- 
domised genomes. In considering the evolution of a self- 
assembling system such as protein quaternary structure 
[7] , it may be that the uniform population of trivial 
phenotypes afforded by our aforementioned zero initial 
conditions is more biologically relevant. 

To compare the two scenarios, we ran simulations of 
the an d ^6,8 systems with random, rather than zero, 
initial conditions. The results (Fig. [TTj ) show a signif- 
icant departure from our results with zero initial con- 
ditions. The difference is particularly pronounced at 
high N, where the diversity provided by a large popu- 
lation of random genomes will lead to very low discov- 
ery times, as space can be explored very quickly from 
this start point before any adaptation takes place. In 
fact, the N = 640 6>2,s system shows a discovery time of 
one, as the proportion of search space corresponding to 
a solution (35 328/2 24 ~ 2.1 x 10" 3 ) is more than 1/N 
(1/640 ~ 1.6 x 10 -3 ), making it likely that at least one 
random genome in the initial population will already be 
a suitable solution. By contrast, this random search ef- 
fect has little impact in the much larger search space of 
the £>6,8 system. 



E. Recombination 

We next set R — 1, modelling sexual reproduction. 
This parameterisation was observed to have little effect 
on the behaviour of r values in the £2,8 an d ^6,8 sys- 
tems with zero initial conditions, leading only to a slight 
increase in adaptation times for given fi. The effect of 
setting R = 1 with random initial conditions was much 
more pronounced. In this case, discovery times were sig- 
nificantly reduced and adaptation times were raised in 
both systems, suggesting that recombination may act to 
increase the 'effective mutation rate' experienced by a 
genome. 

In this picture, recombination may act to decorrelate 
an offspring from both its parents if the genetic diversity 
in the population is high. This effect may be, to first or- 
der, absorbed into an effective mutation rate dependent 
on the diversity in the population. Random initial con- 
ditions ensure that this diversity is high, particularly for 
large N, and hence the steps across genome space caused 
by crossover may be large. This 'genetic drift' acts in 
cohort with the bare mutation rate /i, facilitating rapid 
discovery of solutions on the small £2,8 search space, but 
acting to hinder adaptation at higher fi. 



Many studies of evolution employ random initial con- 
ditions, where the initial population is randomised be- 
fore numerical simulation [26l |50l [51]. While this pic- 
ture is appropriate for the modelling of randomly dis- 
tributed alleles in a population, it is of dubious biological 



F. Elitism 

Optimisation-oriented applications of GAs often em- 
ploy elitism. In a population of N individuals with 
elitism e (where e g [0, 1)), the fittest eN individuals are 
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FIG. 12: (colour online) Adaptation time t (solid lines) and dis- 
covery time tjj (dashed lines) with t = 0.1 for (a) 52,8, (b) £6,8, 
with zero initial conditions. The increase in to with high /x for 
nt = 6 has vanished, and all r values are lower than the e = 
equivalents. 



preserved totally intact from one generation to the next, 
immune to the action of mutation and recombination. In 
this way, the information within the fittest individuals 
— the location of the highest peak thus far discovered — 
is preserved, so that decorrelation from this point pro- 
gresses more slowly and can never be complete. This 
approach is often beneficial for optimisation as it allows 
larger /i values to be used — increasing exploration effi- 
ciency — without loss of information about the current 
best solution. 

The biological relevance of elitism is questionable. The 
problem arises from the immunity of the fittest individu- 
als to mutation (and crossover, in a sexually reproducing 
population). This situation essentially corresponds to a 
number of extremely long-lived individuals which contin- 
ually reproduce through their lifetimes, dying only when 
a fitter solution is found. 

Elitism can have a profound effect on the evolutionary 
dynamics of a model. Fig. 12 shows (/x, r) curves for a 

These ef- 



range of evolutionary scenarios with e = 0.1. 
fects include a general reduction in t values, showing that 
elitism is a useful tool in pure optimisation application of 
GAs. The increase in td with high fj, on 8 is 110 longer 



3 


e 


4 

3 


E 



FIG. 13: Illustration of Di and Ca symmetries in homomeric 
tetramers. 



observed, as elitism retains information from one genera- 
tion to the next, meaning that the search never becomes 
fully random. In experiments with recombination (not 
pictured), elitism also acts to stabilise the population, 
with adaptation observed in e = 0.1 simulations in some 
regimes that struggled to adapt with e = 0. 



V. HOMOMERIC PROTEIN TETRAMERS 

It has been estimated that between 50 and 70% of pro- 
teins form homomeric clusters in vivo |52j . These com- 
plexes are usually symmetrical, with each protein in an 
identical environment. Homomeric tetramers, for exam- 
ple, may display cyclic symmetry (C4) or dihedral sym- 
metry {D2). The C4 geometry involves only one type 
of interaction, whereas the D2 complex involves at least 
two self-complementary interactions. In an important re- 
cent study by Levy et al. [5], it was shown that dihedral 
complexes are over 10 times more abundant than cyclic 
complexes with the same number of subunits. More- 
over, these authors found that the evolutionarily older 
interactions are typically stronger than the more recently 
evolved patches, and that the clusters dissassembled in a 
hierarchical fashion, with the newer (and weaker) bonds 
breaking first. 

The relationship between the strength of the patches 
and their evolutionary history, as well as the observed 
hierarchical dissassembly can be rationalized with simple 
statistical mechanical models [7J. Similarly, the prefer- 
ence for dihedral over cyclic symmetry has been linked 
to the fact that for D2 structures, two pairs of identical 
edges bond (requiring self-complementary interactions or 
homointeractions) , whereas in C4 structures, one pair of 
different edges bond (using non-self-complementary in- 
teractions or heterointeractions) . Statistical models of 
the formation of homointeractions and heterointeractions 
have shown that the former have a wider distribution of 
energies than the latter. It has been suggested that this 
wide distribution makes stable low-energy bonds easier 
to evolve using homointeractions than heterointeractions 
which may result in a biological preference for D2 struc- 
tures [53rl55j . Another reason for the preference for D2 
may be that evolution does not need to proceed to a 
tetramer structure in a single step, but can go through a 
dimeric intermediate. 
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Our simple polyomino model cannot be used in its 
current form to study the strength of patches, and 
by extension, the hierarchical assembly /disassembly. 
However, it can be used to investigate the effect of 
homo/heterointeractions and evolutionary intermediates 
on the evolutionary preference for D2 over C4. In or- 
der to model this system, we must generalise our model 
to allow tetrameric structures to form in both symme- 
try configurations, as shown in Fig. 



13 To do this, we 



allow building block tiles to 'flip', so that, for example, 
tiles {1,2,3,4} and {1,4,3,2} are equivalent. The sides 
of building blocks now correspond to free, rather than 
fixed, necklaces [3T]. This condition reflects the fact that 
homointeraction interfaces require a rotation by ir radi- 
ans with respect to each other to form a bond. 

We first investigate the case where heterointeractions 
are equally easy to evolve as homointeractions. To 
achieve this, we choose a new interaction matrix such 
that the bonding pairs are: 3(4 3,4H4,2f>6, with 
all other colours neutral. This setup was chosen so that, 
given zero initial conditions, the formation of two self- 
interacting edges involves the same number of mutations 
as the formation of a non-self-interacting bonding pair. 
Specifically, the discovery of colours 3 and 4 (Oil and 
100) or 2 and 6 (010 and 110) are equally likely, each 
requiring three beneficial mutations. We label this new 
search space 5'i8, with a characteristic number of self- 
interactions n S i = 2. We note that, given that n t — 1 for 
this system, there is no distinction between the SFN and 
SGN cases mentioned in Section III Al 

In a similar manner to that used for the £2,8 system 



in Section III we can evaluate all possible structures in 
this new search space and the possible transitions be- 
tween phenotypes (see Fig. 14). There are 4096 differ- 
ent possible genotypes, which arc distributed among the 
possible phenotypes as shown in Table [TJ A completely 
random search would thus display a D2 structure fre- 
quency of 0.68. While the interactions are chosen so that 
the minimal number of mutations required to reach a D2 
structure from zero initial conditions is the same as that 
required to reach a C4 structure, the redundancy avail- 
able to D2 genomes (which may contain, for example, one 
unpaired heterointeraction in addition to their homoin- 
teractions) gives -D4 a higher search space volume than 
that of the less redundant C4 structures. We can also 
map out all the pathways between different phenotypic 
states, as shown in Fig. 
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To study the dynamic effects of the structure of search 
space, we simulated a population of 10 5 random walk- 
ers in genome space. Each walker started from zero ini- 
tial conditions, and then took mutational steps until a 
genome encoding one of the two tetrameric states was 
reached. A mutational step involved an application of 
the mutation operator from a GA, rather than enforcing 
exactly one mutation per step. Walks were terminated 
and ignored if they reached the UND state (something 
that mirrors what might happen in nature where this 
usually would be lethal for the organism). 
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FIG. 14: (colour online) Transition probabilities for 5' 1,8. (a) 
Number of self-interacting colours n s i = 2 and (b) n S i = 4. (i) 
Transition probabilities between phenotype x and phenotype 
y. (ii) Transition probabilities represented in a network be- 
tween phenotypes. The edge widths are proportional to their 
probability. 
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TABLE I: Number of genomes in the S' ls search spaces that 
encode different structures. 



A similar random walker analysis is possible in pheno- 
type space, on the network in Fig. 
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ii). Here, each 
random walker occupies a node in the network, and may, 
at each time-step, undergo a transition between nodes 
according to the weight of the connecting edge. A popu- 
lation of walkers was initialised at the monomer node and 
allowed to walk, with UND encounters being terminated 
and ignored. The results of both these walker simulations 
are shown in Table |TTJ 

To test the effect of a dimer intermediate on the prob- 
ability of obtaining a D2 or a C4 structure, we also used 
a GA to run an evolutionary simulation. We employed 
two different fitness functions, representing two situa- 
tions: dimers possessing either no fitness advantage or a 
large fitness advantage over monomers. As the only pos- 



13 



SSP GW PW GA, FF A GA, FF B 

n 3i = 2 0.68 0.58 0.48 0.55 0.67 
n ai = 3 0.88 0.75 0.67 0.69 0.80 
n 3i = 4 0.95 0.86 0.75 0.96 0.98 

TABLE II: Proportion of random walker and GA simulations 
on S[ a that result in a D2 structure being discovered before 
a C4 structure. Columns are Search Space Proportion (SSP), 
defined as the number of genomes encoding D2 structures 
divided by the total number of tetramer genomes, Genotype 
Walker (GW), Phenotype Walker (PW), and GAs, with FF 
denoting fitness function, as described in the text. n si is the 
number of self-interacting colours in the rule set. GAs were 
run with N = 80, fiL = 0.5. 



sible phenotypes in this landscape are UND and s = 1, 
s = 2, s = 4, we represent a fitness function with the 
values awarded to these four cases respectively. Fit- 
ness function A gives no advantage to dimer formation: 
F(UND) = 0,F(1) = 0.1, F(2) = 0.1, F(4) = 1. Fitness 
function B gives a large fitness advantage to dimer for- 
mation: F(UND) = 0,F(1) = 0.1, F{2) = 0.9, F(4) = 1. 
We ran 10 4 simple GAs for each case, with N = 80 and 
[ih = 0.5, and measured the proportion of times a run 
discovered (rather than adapted to) either a C4 or a D2 
phenotype. Table [Tl*| shows the results of simulations with 
these fitness functions. 

We then introduced an evolutionary bias towards ho- 
mointeractions by allowing more colours to self-interact. 
To this end, we first include 1 •<-> 1 (giving n S i = 3) 
and then 5 O 5 bonds (giving n S i = 4). This addition 
of homointeractions changes the evolutionary landscape 
dramatically (see Fig. 14 b)). The distribution of pheno- 
types is shown in Table|IJ The number of UND genotypes 
is observed to increase with n S i, due to the greater num- 
ber of genomes that encode extended, unbound struc- 
tures in systems with large numbers of self-interactions. 

Table [TT] shows a comparison of the D2 : C4 ratio ex- 
pected from search space structure with results for walker 
and GA simulations on these systems. A number of 
trends can be observed in Table [TT] Although the inter- 
actions are chosen so that the minimum number of mu- 
tations required to reach a D 2 structure from zero initial 
conditions is the same as that for a C4 structure, D2 
structures appear more frequently, which is commensu- 
rate with the fact that they occupy a larger proportion 
of search space. 

However, the proportion of runs that first discover D2 
structures is lower than expected from the search space 
structure for genotype walkers, and lower still for pheno- 
type walkers. The slightly lower proportion for genotype 
walkers is due to the starting point of the simulations: 
monomers, of all possible phenotypes, display the highest 
transition probability to C4 structures, so C4 discovery 
is more likely from the zero initial conditions we employ 
(encoding a monomer) than from a random start point. 

The significantly lower D2 proportion from phenotype 



walkers is due to the shorter length of the monomer 
— >• tetramer pathways, which requires only one tran- 
sition, whereas monomer — » dimer — » D 2 requires 
two. In this case, the phenotype representation has 
masked the genetic detail whereby the minimal num- 
ber of steps required to reach D2 and C4 structures 
from zero initial conditions are identical. The steps in 
the minimal monomer — > C4 pathway involve one neu- 
tral monomer step (0000 —> 0002) and one phenotype- 
changing monomer — > C4 step (0002 — > 0062), whereas 
both steps to reach a D 2 structure are phenotype- 
changing (0000 -> 0003 0043). The observed dif- 
ference is an illustration of the influence of a complex 
genotype-phenotype map. In this case, information is 
lost when mutational steps across a neutral network are 
disregarded. 

The proportion of GA runs with fitness function A that 
identify a D 2 structure is lower than expected in com- 
parison to the genotype walkers for n s j = 2 and n S i — 3. 
This difference arises from the different amounts of time 
required for a GA to identify D 2 and C4 structures. For 
n s i = 2 and 3, it is observed that the mean discovery 
time for C4 structures is lower than the mean discovery 
time for D2 structures. A GA reports the structure it 
first discovers, whereas a set of genotype walkers reports 
the proportion of structures discovered regardless of the 
relative time taken to reach these structures. The lower 
C4 mean discovery time for low n S i GAs therefore results 
in more C4 structures being reported than in the geno- 
type walkers. For n S i = 4, the mean discovery time for 
D2 structures is lower than that for C4 in GAs, reflected 
in the higher observation of D2 structures in these GA 
simulations. Note that if the GWs were run in parallel 
sets, and the set was stopped at the first discovery of a 
tetramer, this would also favour C4 for n S i = 2,3 and D2 
for n sl = 4. 

Another effect that acts to change the expected D2 ■ C4 
ratio arises from UND structures. In a GA, genomes en- 
coding UND structures will be replaced (due to their zero 
fitness) by a copy of another genome chosen by selection. 
This replacement genome will be either a monomer or a 
dimer, according to the current state of the GA popula- 
tion. As n s i increases, or if fitness function B is used, the 
population becomes more likely to contain dimers, due 
respectively to their increased presence in search space 
and their increased fitness. If UND genotypes are re- 
placed by monomers, C4 discovery will be more likely 
(the case at low n si ). If they are replaced by dimers, D 2 
discovery will be more likely (the case at high n S i). 

Another noticeable result is that conferring a fitness 
advantage to dimers increases the proportion of D2 struc- 
tures discovered in GAs. This increase is due to selection 
favouring dimers in the evolving population, from which 
situation the dimer — > D 2 transition is most likely. 

The above GA results concern the discovery of 
tetramcrs rather than adaptation of the population to 
tetramers. When adaptation was considered, the n si = 2 
and 3 trends remained very similar. The n S i = 4 sys- 
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tern became 100% dominated by D2 tetramers, as the 
genomes encoding C4 structures in this system were in- 
dividual and isolated. In other words, they exhibit low 
phenotypic robustness and adaptation proved impossible 
with such small neutral network sizes. 

We note that the coarse-grained nature of our model 
greatly simplifies the description of protein surfaces. In 
proteins, interacting sites consist of multiple amino acid 
residues, rather than a single colour type as we employ. 
Point mutations in reality will normally alter not more 
than one constituent amino acid of a bonding site, rather 
than entirely changing the bonding characteristics of an 
interaction site. In addition, the spatial structure of 
protein complexes is vastly more complicated than the 
simple 2D tile geometry we employ here. However, this 
simple system nonetheless displays interesting dynamic 
behaviour. We show that favouring homo-interactions in 
the search space, and favouring dimers in the fitness func- 
tion, can both significantly enhance the proportion of D2 
tetramers over C4 tetramers. By performing a complete 
enumeration of the the fitness landscape, we can also 
uncover some subtle questions related to the underlying 
structure of the landscape. For example, considering only 
the phenotype structure can mask important genotypic 
structure that influences the evolutionary dynamics. 

VI. CONCLUSIONS 

We have studied the evolutionary dynamics of self- 
assembling polyominoes. We focussed on deterministic 
self-assembly - where a given rule set always leads to the 
same polyomino structure - because an analogy can be 
made with monodisperse self-assembly seen in nature, for 
example in protein quaternary structure. 

Although our model is simple enough to be easily 
tractable with modest computational resources, it ex- 
hibits rich evolutionary behaviour that is linked to its 
non-trivial genotype-phenotype mapping. The evolution- 
ary dynamics can be viewed as a search performed by 
a population of individuals on a complex fitness land- 
scape. An advantage of the polyomino system is that in 
some cases this landscape can be fully enumerated and 
classified in terms of adjacent structures and the tran- 
sitions between them. Such information helps explain 
some of the detailed behaviour observed in GA simula- 



tions. Properties like robustness and evolvability [JS] can 
easily be extracted from the fully enumerated landscapes. 

We also investigated the effect of changing the muta- 
tion rate, the population size, and the size of the search 
space on adaptation and discovery times for the evolution 
of certain classes of polyominoes. We find that there is 
an optimal, intermediate mutation rate value for adapta- 
tion. For smaller fj,L the system takes longer to discover 
the desired phenotypes, whereas for larger [iL the mu- 
tational entropy prevents it from adapting to the right 
phenotype. 

For smaller spaces and larger populations the discovery 
time keeps decreasing with increasing mutation rate, but 
for larger spaces, there is also an optimal mutation rate 
for the discovery time. This effect can be cast into the 
language of exploration and exploitation [56-58 ]. For low 
fiL, the system can only take small steps across the search 
space, leading to confinement of the gene pool around 
fitness optima |27j , low diversity, and slow exploration of 
surrounding space. At high /j,L, the system de-correlates 
very quickly, reducing its ability to exploit beneficial mu- 
tants through further small changes, raising diversity to 
almost the level expected for a randomised population. 
The search's hill-climbing ability is decreased as large 
steps randomise the gene pool very quickly. 

The modelling of evolutionary processes with genetic 
algorithms is complicated by the fact that the number of 
parameters that can be varied is very large. One advan- 
tage of our polyomino system is that the effects of varying 
the GA parameters can be easily quantified. We studied 
some popular parameter choices, and argue that, for ex- 
ample, the use of random initial conditions or elitism may 
not be the most biologically relevant way to parameterise 
a genetic algorithm. 

Finally, we studied the evolution of polyomino 
tetramers, inspired by recent work on the structure and 
evolution of homomeric protein tetramers [51 [7j . In na- 
ture there is a strong preference of D2 over C4 symme- 
tries, and we show that both an increase in the proba- 
bility of homointeractions as well as a fitness advantage 
of dimeric intermediates can strongly favour the forma- 
tion of D2 symmetry. Our simplified model shows that 
the outcome of evolutionary dynamics is affected by the 
topology of the search space, including emergent proper- 
ties like phenotypic robustness. 
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